﻿using System;

namespace Atomic.Libraries.Mathematics.Integration
{
	public sealed class DoubleExponentialOscillatory
	{
		/// <summary>
		/// Computes the integral $\int_0^\infty f(x)\Math.Sin\omega xdx$ using the double exponential formula for Fourier transforms.
		/// </summary>
		public static double Integrate(Func<double, double> f, double omega, double h, int n1, int n2)
		{
			if (n1 < 0 || n2 < 0)
			{
				throw new ArgumentException("Invalid number of mesh points.");
			}

			double beta = 0.25;
			double alpha = beta / Math.Sqrt(1.0 + Math.Log(1.0 + Math.PI / h) / (4.0 * h));
			//double alpha = beta / Math.Sqrt(1.0 + Math.Log(1.0 + Math.PI / (omega * h)) / (4.0 * omega * h));
			//alpha = beta / Math.Sqrt(1 + Math.Log((1 + Math.PI / (omega * h)) / (4 * omega * h)));

			double s = 0.0;
			for (int i = -n1; i <= n2; i++)
			{
				double t = i * h;

				/*double phi_x = t / (1.0 - Math.Exp(-2 * t - alpha * (1 - Math.Exp(-t)) - beta * (Math.Exp(t) - 1)));
				// Compute the transformation and its derivative.
				double e1 = Math.Exp(-t);
				double e2 = Math.Exp(t);
				double a = Math.Exp(-2.0 * t - alpha * (1.0 - e1) - beta * (e2 - 1.0));
				double phi, phi1;
				double phi1_num = (Phi(alpha, beta, t + 1e-8) - Phi(alpha, beta, t - 1e-8)) / (2 * 1e-8);
				if (i != 0)
				{
					phi = t / (1.0 - a);
					phi1 = (1.0 - a - t * a * (2.0 + alpha * e1 + beta * e2)) / ((1.0 - a) * (1.0 - a));
				}
				else
				{
					// Take the limit using l'Hôpital's rule.
					phi = 1.0 / (2.0 + alpha + beta);
					phi1 = (1.0 - a - t * a * (2.0 + alpha * e1 + beta * e2)) / ((1.0 - a) * (1.0 - a));
					phi1 = phi1_num;
				}
				if (Math.Abs(phi1 / phi1_num - 1.0) > 1e-4)
				{
					//Console.WriteLine("*");
				}
				if (Math.Abs(phi / phi_x - 1.0) > 1e-14)
				{
					//Console.WriteLine("*");
				}*/

				double phi, phi1 = 0;
				Phi(alpha, beta, t, out phi, ref phi1);
				s += f(Math.PI / (omega * h) * phi) * Math.Sin(Math.PI / h * phi) * phi1;



				//if (double.IsNaN(f(Math.PI / (omega * h) * phi) * Math.Sin(Math.PI / h * phi) * phi1))
				//{
				//    Console.WriteLine("**");
				//}
			}

			return Math.PI / omega * s;
		}

		public static double Integrate(Func<double, double> f, double omega, double eps)
		{
			// http://reference.wolfram.com/mathematica/tutorial/NIntegrateIntegrationStrategies.html#48895312
			throw new NotImplementedException();
		}

		public static void Test()
		{
			//Console.WriteLine(Integrate(x => Math.Log(x), 1.0, 0.05, 100));
			//Console.WriteLine(Integrate(x => x != 0.0 ? 1.0 / Math.Sqrt(x) : 0.0, 1.0, 0.00075, 5000 * 9400, 2000 * 6900));
			//Console.WriteLine(Integrate(x => x != 0.0 ? 1.0 / Math.Sqrt(x) : 0.0, 1.0, 0.075, 94, 69));
			//Console.WriteLine(Math.Sqrt(Math.PI / 2));
			new DoubleExponentialOscillatory().main();
		}

		public static double Phi(double alpha, double beta, double t)
		{
			double e1 = Math.Exp(-t);
			double e2 = Math.Exp(t);
			double a = Math.Exp(-2.0 * t - alpha * (1.0 - e1) - beta * (e2 - 1.0));
			return t / (1.0 - a);
		}

		public static void Phi(double alpha, double beta, double t, out double phi, ref double phi1)
		{
			double e1 = Math.Exp(-t);
			double e2 = Math.Exp(t);
			double a = Math.Exp(-2.0 * t - alpha * (1.0 - e1) - beta * (e2 - 1.0));

			if (t != 0.0)
			{
				phi = t / (1.0 - a);
			}
			else
			{
				// Take the limit using l'Hôpital's rule.
				phi = 1.0 / (2.0 + alpha + beta);
			}

			if (phi1 == 0)
			{
				phi1 = (1.0 - a - t * a * (2.0 + alpha * e1 + beta * e2)) / ((1.0 - a) * (1.0 - a));

				if (double.IsNaN(phi1))
				{
					double h = 1e-8;
					double f, f1 = 1, g, g1 = 1;
					Phi(alpha, beta, t - h, out f, ref f1);
					Phi(alpha, beta, t + h, out g, ref g1);
					phi1 = (g - f) / (2.0 * h);
				}

			}

		}

		// http://www.kurims.kyoto-u.ac.jp/~ooura/profile.html
		// http://www.kurims.kyoto-u.ac.jp/~ooura/papers/deft_new.c

		double M_PI = 3.14159265358979323846;
		double M_EULER = 0.57721566490153286060;

		void phi2init(double epsm, double epsp, double h, double o0, out int nm, out int np)
		{
			double M, a, b, t;

			M = M_PI / (o0 * h);
			b = 0.25;
			a = b / Math.Sqrt(1 + M * Math.Log(1 + M) / (4 * M_PI));
			t = Math.Log(-Math.Log(epsp / M) / b + 1);
			t = Math.Log((2 * t - Math.Log(epsp / M / t)) / b + 1);
			np = 1 + (int)(Math.Log((2 * t - Math.Log(epsp / M / t)) / b + 1) / h);
			/* |M*(*np*h - phi(*np*h))| = O(epsp) */
			t = Math.Log(-Math.Log(epsm / M) / a + 1);
			t = Math.Log((2 * t - Math.Log(epsm / M / t)) / a + 1);
			nm = -1 - (int)(Math.Log((2 * t - Math.Log(epsm / M / t)) / a + 1) / h);
			/* |M*phi(*nm*h)| = O(epsm) */
		}

		void phi2(int nm, int np, double h, double o0, NonZeroArray aw)
		{
			int n;
			double M, a, b, t, em, ep, u;

			M = M_PI / (o0 * h);
			b = 0.25;
			a = b / Math.Sqrt(1 + M * Math.Log(1 + M) / (4 * M_PI));
			for (n = nm; n <= np; n++)
			{
				t = h * n;
				if (n != 0)
				{
					em = a * (1 - Math.Exp(-t));
					ep = b * (Math.Exp(t) - 1);
					u = 1 / (1 - Math.Exp(-2 * t - em - ep));
					aw[2 * n] = M * t * u; /* M * phi(n*h) */
					aw[2 * n + 1] = M * h * (1 - t * (2 + a - em + b + ep) * (u - 1)) * u; /* M*h * phi'(n*h) */
				}
				else
				{
					u = 1 / (2 + a + b);
					aw[2 * n] = M * u; /* M * phi(0) */
					aw[2 * n + 1] = M * h * (1 + (a - b) * u * u) / 2; /* M*h * phi'(0) */
				}
			}
		}

		void deft2init(int nm, int np, double o0, int komega, NonZeroArray aw, NonZeroArray wk, Func<double, double> f)
		{
			int n;
			double ostep, x, xi, y, sr, si, wr, wi;

			ostep = o0 / komega;
			for (n = nm; n <= np; n++)
			{
				x = aw[2 * n];
				sr = Math.Cos(ostep * x);
				si = Math.Sin(ostep * x);
				wk[5 * n] = 2 * si;
				xi = (x * o0 - M_PI * n) * 0.5;
				wi = Math.Cos(xi);
				wr = Math.Sin(xi);
				y = aw[2 * n + 1] * 2 * wr * f(x);
				wr *= y;
				wi *= y;
				wk[5 * n + 1] = wr;
				wk[5 * n + 2] = wi;
				wk[5 * n + 3] = wr * sr - wi * si;
				wk[5 * n + 4] = wr * si + wi * sr;
			}
		}

		void deft2(int nm, int np, int komega, NonZeroArray wk, double[] F)
		{
			int k, n;
			double s1r, s1i, s2r, s2i, ws, w1r, w1i, w2r, w2i;

			for (k = 0; k < komega; k++)
			{
				s1r = 0;
				s1i = 0;
				s2r = 0;
				s2i = 0;
				for (n = nm; n <= np; n++)
				{
					ws = wk[5 * n];
					/*
					s1r += wk[5 * n + 1];
					s1i += wk[5 * n + 2];
					s2r += wk[5 * n + 3];
					s2i += wk[5 * n + 4];
					wk[5 * n + 1] -= ws * wk[5 * n + 4];
					wk[5 * n + 2] += ws * wk[5 * n + 3];
					wk[5 * n + 3] -= ws * wk[5 * n + 2];
					wk[5 * n + 4] += ws * wk[5 * n + 1];
					*/
					w1r = wk[5 * n + 1];
					w1i = wk[5 * n + 2];
					s1r += w1r;
					s1i += w1i;
					w2r = wk[5 * n + 3];
					w2i = wk[5 * n + 4];
					s2r += w2r;
					s2i += w2i;
					w1r -= ws * w2i;
					w1i += ws * w2r;
					wk[5 * n + 1] = w1r;
					wk[5 * n + 2] = w1i;
					wk[5 * n + 3] = w2r - ws * w1i;
					wk[5 * n + 4] = w2i + ws * w1r;
				}
				F[4 * k] = s1r;
				F[4 * k + 1] = s1i;
				F[4 * k + 2] = s2r;
				F[4 * k + 3] = s2i;
			}
		}

		double f1(double x)
		{
			return Math.Log(x) / Math.Sqrt(x);
		}

		double F1r(double o)
		{
			return Math.Sqrt(M_PI / 2 / o) * (Math.Log(4 * o) + M_EULER + M_PI / 2);
		}

		double F1i(double o)
		{
			return Math.Sqrt(M_PI / 2 / o) * (Math.Log(4 * o) + M_EULER - M_PI / 2);
		}

		void main()
		{
			int np, nm, komega, k;
			double h, o0, o, er, ei;
			NonZeroArray aw = new NonZeroArray(2000, -1000);
			NonZeroArray wk = new NonZeroArray(5000, -2500);
			double[] F = new double[1000];

			h = 0.075;
			o0 = 1;
			komega = 128;
			phi2init(1.0e-24, 1.0e-12, h, o0, out nm, out np);
			phi2(nm, np, h, o0, aw);
			deft2init(nm, np, o0, komega, aw, wk, f1);
			deft2(nm, np, komega, wk, F);
			for (k = 1; k < 2 * komega; k++)
			{
				o = o0 * k / komega;
				er = F[2 * k] + F1r(o);
				ei = F[2 * k + 1] + F1i(o);
				Console.WriteLine("omega= {0}, error= {1}", o, Math.Sqrt(er * er + ei * ei));
			}
			Console.WriteLine("N_-= {0}, N_+= {1}", nm, np);
		}

		private class NonZeroArray
		{
			private double[] values;
			private int offset;

			public NonZeroArray(int length, int offset)
			{
				this.offset = offset;

				values = new double[length];
			}

			public double this[int i]
			{
				get
				{
					return values[i - offset];
				}
				set
				{
					values[i - offset] = value;
				}
			}
		}
	}
}
